######################################################################
#
#Michael Poznansky and Matt K. Scroggs
#
#"Ballots and Blackmail: Coercive Bargaining and the Democratic Peace"
#International Studies Quarterly
#                                                                           
######################################################################


#Purpose
#This R file walks through the process of creating the standard errors for the cluster robust sandwich models.
#
# https://dataverse.harvard.edu/dataverse/mkscroggs
#
# Version 1.0
# Last updated: March 29, 2016
library(foreign)
install.packages("sandwich")
install.packages("lmtest")
library(sandwich)
library(Matrix)
rm(list=ls())

# Cluster robust variance estimation function
robust.se.nodfc <- function(model, cluster){
	require(sandwich)
	require(lmtest)
	M <- length(unique(cluster))
	N <- length(cluster)
	K <- model$rank
	dfc <- 1
	uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum))
	rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
	rcse.se <- coeftest(model, rcse.cov)
	return(list(rcse.cov, rcse.se))
}

##############
#PAPER MODELS#
##############

#INPUT VARIABLES FOR TABLE 3, MODEL 6 (polityjoint)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/lowdpnd_sandwich.dta")
regSpec <- paste("mct~", 
				paste("polyprod","lowdpnd", "cincratio","alliance1","absidealdiff", "conttype2",
					"time","time_sq","time_cu",
				sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
	iUp <- index[i]
	clusUp <- apply(dyad.mat, 
					1,
					function(x)as.numeric(iUp %in% x))
	clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
	if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
	if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
	cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polyprod lowdpnd cincratio alliance1 absidealdiff conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#INPUT VARIABLES FOR TABLE 4, MODEL 6 (polity_low)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/lowdpnd_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_low","lowdpnd", "cincratio","alliance1","absidealdiff", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_low lowdpnd cincratio alliance1 absidealdiff conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#INPUT VARIABLES FOR TABLE 5, MODEL 6 (polity_both)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/lowdpnd_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_both","lowdpnd", "cincratio","alliance1","absidealdiff", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_both lowdpnd cincratio alliance1 absidealdiff conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#################
#APPENDIX MODELS#
#################

#UN AFFINITY (polityjoint)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/UNaffinity_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polyprod","lowdpnd", "cincratio","alliance1","UN_affinity2", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polyprod lowdpnd cincratio alliance1 UN_affinity2 conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#UN AFFINITY (polity_low)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/UNaffinity_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_low","lowdpnd", "cincratio","alliance1","UN_affinity2", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_low lowdpnd cincratio alliance1 UN_affinity2 conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#UN AFFINITY (polity_both)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/UNaffinity_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_both","lowdpnd", "cincratio","alliance1","UN_affinity2", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_both lowdpnd cincratio alliance1 UN_affinity2 conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#LOW CIE (polityjoint)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/lowCIE_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polyprod","lowdpnd", "lowCIE", "cincratio","alliance1","absidealdiff", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polyprod lowdpnd lowCIE cincratio alliance1 absidealdiff conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#LOW CIE (polity_low)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/lowCIE_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_low","lowdpnd", "lowCIE", "cincratio","alliance1","absidealdiff", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_low lowdpnd lowCIE cincratio alliance1 absidealdiff conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#LOW CIE (polity_both)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/lowCIE_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_both","lowdpnd", "lowCIE", "cincratio","alliance1","absidealdiff", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_both lowdpnd lowCIE cincratio alliance1 absidealdiff conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#UN AFFINITY & LOW CIE (polityjoint)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/UNAffinityCIE_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polyprod","lowdpnd", "lowCIE", "cincratio","alliance1","UN_affinity2", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polyprod lowdpnd lowCIE cincratio alliance1 UN_affinity2 conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#UN AFFINITY & LOW CIE (polity_low)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/UNAffinityCIE_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_low","lowdpnd", "lowCIE", "cincratio","alliance1","UN_affinity2", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_low lowdpnd lowCIE cincratio alliance1 UN_affinity2 conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64

#UN AFFINITY & LOW CIE (polity_both)

triangle <- read.dta("~/Dropbox/Interdependence--CompellentThreats/IO-ISQ/ISQ Note/Accepted Version/Datasets/UNAffinityCIE_sandwich.dta")
regSpec <- paste("mct~", 
                 paste("polity_both","lowdpnd", "lowCIE", "cincratio","alliance1","UN_affinity2", "conttype2",
                       "time","time_sq","time_cu",
                       sep="+"), sep="")

regForm <- as.formula(regSpec)
fit <- glm(regForm, data=triangle, x=T, family="binomial")

index.sc <- unique(c(triangle$statea, triangle$stateb))
index <- index.sc[order(index.sc)]
dyad.mat <- triangle[,c("statea","stateb")]

# Dyadic cluster robust via multiway decomposition
for(i in 1:length(index)){
  iUp <- index[i]
  clusUp <- apply(dyad.mat, 
                  1,
                  function(x)as.numeric(iUp %in% x))
  clusIndexUp <- clusUp*-99 + (1-clusUp)*1:nrow(dyad.mat)
  if(i==1){dcrUp <- robust.se.nodfc(fit, clusIndexUp)[[1]]}
  if(i>1){dcrUp <- dcrUp + robust.se.nodfc(fit, clusIndexUp)[[1]]}	
  cat(paste(iUp, " ")); flush.console()
}
# substract naive CR:
dcrUp2 <- dcrUp - robust.se.nodfc(fit, triangle$dyadid)[[1]]
# substract HR:
Vhat <- dcrUp2 - (length(index)-2)*diag(vcovHC(fit,type="HC0"))
sqrt(diag(Vhat))

# (Intercept) polity_both lowdpnd lowCIE cincratio alliance1 UN_affinity2 conttype2 time time_sq time_cu 

abs(coef(fit)/sqrt(diag(Vhat))) < 1.64